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Abstract 

The two-dimensional kinetic Ising model, when exposed to an oscillating applied magnetic field, 
has been shown to exhibit a nonequilibrium, second-order dynamic phase transition (DPT), whose 
order parameter Q is the period-averaged magnetization. It has been established that this DPT falls 
in the same universality class as the equilibrium phase transition in the two-dimensional Ising model 
in zero applied field. Here we study for the first time the scaling of the dynamic order parameter 
with respect to a nonzero, period-averaged, magnetic 'bias' field, Hf,, for a DPT produced by a 
square-wave applied field. We find evidence that the scaling exponent, 5d, of Hb at the critical 
period of the DPT is equal to the exponent for the critical isotherm, S e , in the equilibrium Ising 
model. This implies that is a significant component of the field conjugate to Q. A finite- 
size scaling analysis of the dynamic order parameter above the critical period provides further 
support for this result. We also demonstrate numerically that, for a range of periods and values 
of Hb in the critical region, a fluctuation-dissipation relation (FDR), with an effective temperature 
T e fj (T, P, Hq) depending on the period, and possibly the temperature and field amplitude, holds 
for the variables Q and Hb- This FDR justifies the use of the scaled variance of Q as a proxy for 
the nonequilibrium susceptibility, d(Q) /dHb, in the critical region. 

PACS numbers: 05.70.Ln, 64.60.Ht, 89.75.Da, 75.70.Cn 
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I. INTRODUCTION 



The dynamic phase transition (DPT) in a ferromagnetic system below its critical temper- 
ature was first observed in numerical solutions of a mean-field model exposed to an oscillating 
mag net 1C aeld QQ. It was then studied tothe, both in mean-field mo d* fiflflB and 

nnQndn 

in kinetic Monte Carlo (KMC) simulations [5|, [6j, [j], [8|, |9j, [10J. A review of this early work can 



be found in Ref . [TlJ . More recently, the study of the DPT has expanded to include varying 
(and often more physical) model geometries. These include mean-field studies of domain- 



wall motion in an anisotropic XY model in one dimension 
of a three-dimensional Ising system 



15 



12 



13 



14j, KMC simulations 



Heisenberg system in an off-axial field 



, and KMC simulations of a uniaxially anisotropic 



161 ]. an el 



with the effect of a thin-film surface energy 
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iptically polarized applied field [l7j, and 
19, [2^]. The phenomenon has also been 



observed in simulations of CO oxidation under oscillating CO pressure 



21 



Further 



simulation studies of the DPT in the two-dimensional kinetic Isin g m odel have appeared 



23, 



24. 



25. 2i 



1Q 



as well as analytical studies of the DPT 



28, 29 



30. 



3l|. 



Here, we concentrate on the DPT in the two-dimensional kinetic Ising model. It was ob- 
served in simulations of this model that there exists a singularity at a critical period of the 
applied oscillating field [s,[l^], and that the critical exponents (3 and 7 (and, with less accu- 
racy, v) are consistent with the universality class of the equilibrium two-dimensional Ising 
transition in zero field [23]] (f3 = 1/8,7 = 7/4, and v = 1). In those studies, the techniques 
of finite-size scaling were extended to the study of the dynamic order parameter (Q, defined 
in Sec. [TT]) in the non-equilibrium steady state. This provided evidence for a diverging cor- 
relation length at a critical value of the period. In particular, because the field conjugate 
to Q and a fluctuation-dissipation relation were not known, a susceptibility could not be 
measured directly, and the scaled variance X® = L 2 ((Q 2 ) — (Q) 2 ), where L is the linear 
system size, was used as a proxy. An analytical argument, based on the correspondence 
of the two-dimensional kinetic Ising model and the continuous, two-dimensional Ginzburg- 
Landau model at the equilibrium critical point, provided an effective Hamiltonian for the 
non-equilibrium system and confirmed that the DPT is in the Ising universality class 28]. 
These findings are consistent with earlier symmetry arguments that any continous phase 
transition in a stochastic cellular automaton that preserves the Ising up-down symmetry 
should be in the equilibrium Ising universality class 32|, l33 |. 
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Recently, experiments were performed on a [Co(0.4nm)/Pt(0.7nm)] 3 multilayer fi 



strong uniaxial anisotropy 



m with 
35, 391. 



34| , whose equilibrium behavior is known to be Ising-like 
The film was exposed to an oscillating (sawtooth) applied field with varying period, in the 
presence of constant 'bias' magnetic fields of varying strength and sign. (The bias field 
is defined explicitly in Sec. [Ill) The behaviors of the dynamic order parameter and its 
variance, as functions of the applied field period and the bias field, provided strong evidence 
for the presence of the DPT in this experimental system. The observed behavior of the order 
parameter with respect to the bias field supported previous conjectures that the conjugate 
field could include the period- averaged magnetic field as an important component, and 
stimulated the numerical investigations in this paper. 

This paper is organized as follows. In Sec. [Ill we describe the two-dimensional kinetic 
Ising model and our computational methods. In Sec. IIII| we verify directly the scaling 
of the dynamic order parameter with respect to the period-averaged magnetic field at the 
critical period, with scaling exponent 5d ~ S e = 15, in agreement with the equilibrium Ising 
transition. In Sees. II VI and IV l we derive the expected asymptotic scaling functions in a finite- 
size scaling analysis of the dynamic phase transition with non-zero period-averaged bias field, 
and then compare the expected scaling of the dynamic order parameter to our numerical 
results. In Sec. IVIj we present numerical data assessing the applicability of a fluctuation- 
dissipation relation (FDR) to this far-from-equilibrium system. We then in Sec. I VIII compare 
the expected scaling of the susceptibility (of the dynamic order parameter) to our numerical 
data, using the results of Sec.EDto reconcile our findings with previous results on the scaling 
of the fluctuations of the dynamic order parameter. Finally, we present a summary of our 
results in Sec. Villi 

II. COMPUTATIONAL MODEL 

In order to facilitate comparison with previous results, we employ the same model and 



computational method as in Ref . [23j . Specifically, we perform kinetic Monte Carlo (KMC) 
simulations of a two-dimensional periodic square lattice of Ising spins Si, which can take 
only the values Si — ±1. The Hamiltonian of the model is 

n = -j^2s i s j -H(t)Y t s i , (i) 

(ij) » 
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where J > is the ferromagnetic exchange interaction, Y^iUj) runs over a h nearest-neighbor 
pairs, Yli runs over all L 2 lattice sites, and H(t) is an oscillating, spatially uniform applied 
magnetic field. The form of H(t) is here taken as a square wave with amplitude H = 0.3J 
and period P, measured in Monte Carlo steps per spin (MCSS). The square- wave form not 
only allows for more efficient KMC simulation, but also reduces the critical period and the 
finite-size effects 
and sawtooth 



or the DPT [23]. Other symmetric field shapes, such as sinusoidal [9l. llOl] 
34|, yield essentially the same results, but with a larger critical period and 
with stronger finite-size effects. The Glauber single-spin-flip MC algorithm with updates at 
randomly chosen sites is used, in which each attempted spin flip is accepted with probability 

W (St -> -Si) = -. An/m „ (2) 

v ' l + exp(A£/T)' v ' 

where AE is the energy change that would result from acceptance of the spin flip, and T 
is the absolute temperature in energy units (i.e., with Boltzmann's constant set to unity). 
All simulations were performed at T = 0.8T C , where T c = 2.269J is the equilibrium critical 
temperature of the square-lattice Ising ferromagnet in zero applied field 37]. 

The system responds to the oscillating field via the time-dependent magnetization per 
site, 

i? 

™w=7ii>( £ )- (3) 

n 

The dynamic order parameter is defined as the average of m(t) over a given field cycle i |l[ : 



1 



iP 



Q* = -5 m{t)dt. (4) 

J(i-i)p 

We define the bias field, so named because it measures the shift (or 'bias') of the periodic 
field toward either negative or positive field values, as the period-averaged magnetic field, 

i r p 

H b = - H(t)dt. (5) 
" Jo 

This definition applies generally to any periodic magnetic field H(t). In this paper, the 
applied field consists of a square wave with period P superposed with a constant magnetic 
field. Applying (jSJ), since the period-average of the square- wave field is zero, the bias field 
H b in this case is simply equal to the superposed constant magnetic field. 
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III. SCALING WITH RESPECT TO THE BIAS FIELD 



In the two-dimensional equilibrium Ising model, the critical isotherm is given (in the 
thermodynamic limit, i.e., as L —>■ oo) as 

m(T = T c ,H ^0)ocH 1/s % (6) 

where the critical exponent S e = 15 j^]. For finite systems, this relationship breaks down 
when the infinite-system correlation length, ^ (T = T c , H) , which diverges as H — > 0, 
becomes comparable to the linear system size L. The relationship also naturally breaks 
down at larger fields away from the critical region. Therefore, a plot of m vs H for a given 
system size L will follow the power law ([6]) for a range of H near the critical value H = 0, 



39) 



with this range extending to smaller H as L is increased 

We can determine directly whether the non-equilibrium system exhibits a similar rela- 
tionship, 

(Q) (P = P c , H b ^0)cx Hl /S * (7) 

in an analogous way. In Fig. [H we plot (Q) vs Hb at the critical value of the period, P = P c . 
In previous work, the reversal time for the magnetization, following instantaneous reversal 
of the uniform magnetic field H at (H = 0.3 J, T = 0.8TA was found as r = 74.5977 MCSS 
[9|, HDJ ■ The critical scaled half-period for the square waveform was determined to be C 



PJ (2r) = 0.918 ± 0.005 [23j. This yields P c = 136.96 ± 0.75 MCSS, and in our simulations 
and analysis in this paper we use P c = 136.96 MCSS. 

A power-law dependence is indeed seen to hold in Fig. [lj within a range which extends 
to lower values of Hf, as L is increased. We fit the L = 256 data between the points 
labeled A and B in Fig. (U finding a statistically significant fit with power-law exponent 
5d = 14.85 ± 0.18. As including points with Hf, > 0.01J was found to greatly reduce the 
statistical significance of the fit, the value Hb = 0.01J serves as a boundary of the scaling 
region at P = P c . This result is consistent with an exponent 5^ = 5 e = 15, suggesting that 
the bias field Hb, for these parameters and the square waveform, is the dominant component 
of a conjugate field which exhibits the same scaling exponent in the DPT as does the applied 
magnetic field in the equilibrium Ising transition. 
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IV. FINITE-SIZE SCALING ANALYSIS WITH BIAS FIELD 



To provide more complete evidence that iff, is the dominant component of the field 
conjugate to (Q), in the next several sections we will demonstrate data collapse onto a two- 
parameter finite-size scaling function for the system-size dependent quantity (Q)l at points 
(P > P c ,Hb > 0), for lattice sizes L = 90, 128, 180, and (in several cases) L = 256, using 
the critical exponents for the equilibrium Ising system. In this section, we briefly review the 
theory of finite-size scaling as it applies to this system. We then determine the expected 
asymptotic forms of the scaling functions, which are compared in later sections of the paper 
to our computational data. 



The theory of finite-size scaling 



40 



41] states that near a continuous phase transition, 
the singular part of the free-energy density for a (i-dimensional system of linear size L can 
be written as 

f L ^L- d Y ± (\e\L^,HL^), (8) 

where e = (T — T c )/T c , H (in units of k B T) is the field conjugate to the order parameter, v 
is the critical exponent for the correlation length, (3 is the exponent for the order parameter, 
5 is the exponent for the critical isotherm, and Y± are scaling functions above (+) and below 
(— ) the critical point. This yields for the order parameter at finite L 

™ L = §t = L-fi/"^ (\ e \L^, HL^) , (9) 

where the exponent for L in the prefactor is obtained by using the hyperscaling relation 
dv = 2 — a and the exponent equality a = 2 — (3(5 + 1). Further differentiation yields the 
susceptibility, 

XL = ^ = L^G 0± (\e\L^,HL^) } (10) 

where the exponent for L in the prefactor is obtained by using the exponent equality 7 = 
(3(6-1). 

It has previously been shown analytically that the DPT for a sinusoidal applied field, 
which is symmetric under H(t) — > —H(t + P/2) and so which can safely be assumed to have 
H c = 0, has an effective Ginzburg-Landau free-energy density in the same universality class 
as the equilibrium Ising model [28j. It therefore appears reasonable to write corresponding 
scaling functions for the dynamic order parameter (Q) and its associated susceptibility x, 

(Q) L = L-V V T ± {\e\L 1 '^ (HJJ)L^) , (11) 



and 

X L = Vl v G± {H c /J)L ps ' v ) , (12) 

where 9 = (P — P c )/P c , and H c is the (as yet unknown) field conjugate to (Q). In this 
paper we express H c (and Hi) in units of the exchange constant, J, so that the second 
scaling parameter in Eqs. [TT] and [12] is dimensionless. The specific form, H c j J, with which 
H c is assumed to enter the second scaling parameter needs more theoretical investigation, 
and could conceivably change as the theory of the DPT is further developed. However, this 



should not affect our conclusions [42j . Computational results for sinusoidal and square- wave 
fields, which both are symmetric under H(t) — > —H(t+P/2) and so presumably have H c = 0, 
have previously confirmed the scaling behavior with respect to 9. The exponent values were 
determined as 7/1/ = 1.74 ± 0.05, (3/v = 0.126 ± 0.005, and v = 0.95 ± 0.15 [23|, consistent 
with the exact values for the two-dimensional equilibrium Ising model, 7 = 7/4 = 1.75, 
j3 = l/8 = 0.125, and v = 1. 

We now determine the expected asymptotic forms of the scaling functions F + (yi, y 2 ) and 
G+(yi, 2/2), where we emphasize that the + subscript indicates that the scaling functions refer 
to the range P > P c , and where the scaling parameters are y\ = 91}I V and y 2 = (H c / J)L@ S I U . 
yj ^> Uz - We expect \l ~ 9~ J = L^^y^ 1 (independent of y 2 ) and (Q)l = XlH c ~ 9~~ t H c ~ 
L~^' v y^ 1 y2i where 7 = [3(8 — 1) was used to obtain the exponent for L in (Q)l- 
yj < y 7 . We expect that (Q) L ~ H l c /S ~ L^l v y X 2 l& and xl = d(Q) L /dH c ~ L^ v y^' S)/S 
(both independent of yi), where 7 = /?(<5 — 1) was used to obtain the exponent for L in 
Thus, the asymptotic forms of the scaling functions are expected to be 

1 V7 J V2 for Vi ^ V2 
F+(yi,y2) = L 0/l/ (Q) L ~ { (13) 



and 



y\ /5 for yi < y 2 



■ -Ui 7 for yi y 2 

g + {y 1 , V2 )=L-T"xL~< %„ 5)/s f ■ ^ 

y\ for < y 2 



V. COMPARISON OF FIRST SCALING FUNCTION TO COMPUTATIONAL 
RESULTS 

In Fig. 12(a), using the equilibrium values (3 e and u e in calculating J-'+fyi, 7/2) = LP' V (Q) l, 
we present a plot of the scaling function JF + vs jyi for different values of y 2 , for lattice 



S 



sizes L = 90, 128, and 180. Here and for the remainder of the paper, exponents with the 
subscripts 'd' and 'e' refer to the behavior of the nonequilibrium system (with a dynamic 
phase transition) and the equilibrium system, respectively. The scaling function exhibits 
a power-law dependence in the regime y\ ^> y 2 , which is consistent with Eq. (1131) . At 
progressively larger values of the constant y 2 , the power-law scaling can be seen to begin at 
increasing values of y±, as would be expected. A best-fit line to the final five points of the 
L = 180 data at y 2 = 3.39 yields an estimate of the scaling exponent — 7d = —1.76 ± 0.07 
in Eq. (ITBl . This is consistent with the previous results for H c = cited above |23j, and 
it supports the hypothesis that jd = 7e = 7/4 = 1.75. In Fig. EJ^b), we present just the 
data for y 2 = 3.39, including additional data points at y\ = 280 and 477. The data deviate 
from the power-law behavior for L = 90 at y\ > 149, for L = 128 at yi > 280, and for 
L = 180 at yi > 477. This locates the boundary of the scaling regime (for y 2 = 3.39) at 
= yi /L l l v « 2.65 . 

In Fig. [HI again using f3 e and v e in calculating J r + , we plot the scaling function JF + vs 
y 2 at different values of yi, in order to examine the scaling behavior for y\ ^> y 2 . For the 
constant values y± = 43.4, 69.7, and 149, power-law scaling can be observed in the regime 
Hi ^ 2/2- A best-fit line to the y\ = 149 data for the five points from y 2 = 3.39 to 84.6 yields 
a scaling exponent of 1.01 ± 0.01, which is consistent with the value of 1 expected from Eq. 

In order to investigate the scaling of T+ in the asymptotic limit y\ y 2 , we plot in Fig. |4] 
the scaling function J 7 + (yi,y 2 ) vs y 2 at the critical period P = P c (i.e., yi = 0), at lattice 
sizes L = 128, 180, and 256. In the range 20 < y 2 < 50, power-law scaling is observed for all 
three lattice sizes. For y 2 > 50, the data deviate from power-law scaling, with the smallest 
lattice size deviating first, as expected in a finite-size scaling plot. A fit of the L = 256 data 
in the range from y 2 = 8.46 to 84.6 produces a scaling exponent 0.0673 ± 0.0008. Since the 
constant factors LP 5 ' 1 ' and L^l v do not affect the fit of the scaling exponent, this is the same 
exponent found in the fit of (Q) vs Hf, in Fig. [U The reciprocal of this scaling exponent 
is thus Sd = 14.85 ± 0.18, which is consistent with the exponent of the critical isotherm, 
5 e = 15, in the equilibrium Ising model. 

The comparison of the second scaling function, G+(yi,y 2 ), to numerical data is more 
clearly presented after the relationship of the susceptibility \l and the scaled variance X® 
has been examined. Therefore, we present in the next section numerical results on the extent 
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of applicability of an FDR between xl and X^, before turning in Sec. IVHI to the second 
scaling function. 



VI. APPLICABILITY OF A FLUCTUATION-DISSIPATION RELATION 

FDRs, such as the Einstein relation, Green-Kubo relations, etc., hold a central place in 
equilibrium statistical mechanics. This is essentially a consequence of detailed balance and 
the role of the partition function as a moment-generating function, and thus such relations 
cannot be readily extended to nonequilibrium steady states. However, it has recently been 
shown that certain FDRs can be extended to far-from equilibrium steady states by use of an 



effective temperature 43|, |44( . Here we will therefore consider whether the nonequilibrium 



susceptibility and the scaled variance of the dynamic order parameter can be related as 

5 _d(Qh L*{{(?) L -{Q)l) _X Q L 

XL = ~mr = t = t~> (15) 

Orib J-eS J-eS 

with an effective temperature T e g, in a way analogous to the equilibrium FDR, 

_ 9{m) L L 2 ((m 2 ) L - (m) 2 L ) 
XL = ^ jr = f , (16) 

in which T is the temperature. As mentioned in Sees. [J and IIVI this conjecture motivated 
the use in previous work of the scaled variance X® as a proxy for xl in investigating the 
scaling behavior of the nonequilibrium system near its critical period. 

To test the extent to which Eq. (fl5|) holds, we computed values of xl and X® for a range 
of periods from P = 140 to 250 MCSS and a range of bias fields from H b — to an upper 
limit between 0.005J and 0.2J. (The bias field necessary to 'saturate' the nonequilibrium 
system, i.e., to produce values of xl and X® near zero, increases as the period is increased.) 
The computations were perfomed at L = 180. The quantity xl was computed directly as a 
numerical derivative: 

Xl(P, Hb) « ((Q)(P, H b + AH b ) - (Q)(P, H b - AH b )) /2AH b . (17) 

The choice of AH b = 0.1H b was found to produce sufficiently accurate values of the numerical 
derivative across the range of bias fields studied. The results for periods P = 140 through 
190 MCSS are shown in Fig. [51 A linear relationship is seen to exist between xl and X®, for 
each value of P, over a wide range of xl values. At each period, the dependence becomes 
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nonlinear below a certain value of xl, as illustrated for periods P = 150, 170, and 190 MCSS 
in Fig. [6j Since low values of the susceptibility xl correspond to large values of the bias 
field Hb, we interpret this breakdown of linearity as an indication that the FDR in Eq. (TIB]) 
holds only in a scaling regime around the critical point, i.e., for a limited range of Hi around 
H b = 0. 

The relationship between X® and xl at the higher periods, P = 220 and P = 250 MCSS, 
is more complicated, as shown in Fig. [7J At P — 220 MCSS, following the nonlinear regime 
at low xl, there is a linear relationship with slope T eS = (6.27 ± 0.11) J up to xl ~ 13 J -1 , 
followed by a second distinct linear dependence with slope T e g « (4.16 ± 0.29) J above 
Xl — 13 J -1 . At P = 250 MCSS, the initial nonlinear dependence is again present. Then 
the first linear regime has T e g = (6.49 ± 0.07) J up to xl ~ 13 J -1 , and is followed by a 
regime which can be characterized as either linear with very gentle slope (0.27 ± 0.22) J, or 
as an effective 'saturation' of X® past xl = 13 J~ l - 

In Fig. [8] we plot the best-fit slopes from Figs. [5] and U\ which according to Eq. ( fl5l) 
represent estimates of T e s, vs the scaling parameter 9 = (P — P c ) / P c . We have included 
in the plot the slopes of both linear regimes for the values 9 = 0.606 and 0.825 (P = 220 
and 250 MCSS). For 9 below 0.4 (P « 190 MCSS), T eff increases with in a way not 
inconsistent with a linear relationship (with slope 2. 97 J). It may be interesting to note that 
an extrapolation of the linear relationship to 9 = (P = P c ) yields the value T e g = 3.39J, 
which is significantly higher than the critical temperature, T c = 2.2619 J, of the equilibrium 
Ising system. However, one should not put too much emphasis on the numerical values of 
T e s 45] , as they could easily be changed. For instance, if Hb is only proportional to the full 
conjugate field H c , with a proportionality constant different from unity, this would trivially 
change T e g in Eq. ([15]) . The important result, which we have demonstrated to hold in the 
critical region, is the linear relationship between X® and xl- 

We can thus characterize the extent of applicability of an FDR to the DPT above the 
critical period as follows. For 9 < 0.4, an FDR holds outside of a small nonlinear regime at 
low xl (high H^), with an effective temperature T e ff which increases approximately linearly 
with 9. For 9 above 0.4, two linear relationships appear to exist between X® and xl i n 
separate regimes, making it impossible to define a unique T e g at a given value of 9. An 
understanding of the nonlinear regime, which is present at low xl for all periods examined, 
as well as of the complication of the FDR above 9 = 0.4, would be highly desirable. We hope 



11 



that these numerical results can stimulate the development of, as well as test the accuracy 
of, a theoretical description of the non-equilibrium steady states produced in the presence 
of non-zero for this DPT. 

VII. COMPARISON OF SECOND SCALING FUNCTION TO COMPUTA- 
TIONAL RESULTS 

We now test the asymptotic scaling forms for Q + in Eq. (|14|) . First, we note that in 
performing least-squares fits, one normally requires the goodness-of-fit parameter q, i.e., the 
probability that (assuming the fit relationship were true) random error alone could produce 
the observed data, to be greater than 10~ 3 to consider the fit reasonable. Within this section, 
however, and in the captions to Figs. [9] through [Til h wm De useful for descriptive purposes 
to refer to scaling exponents resulting from attempts at least-squares fits with q values below 
this acceptable range. We will refer to the results of such unsuccessful fitting attempts as 
'nominal' scaling exponents, and for clarity will report the value of the parameter q for each 
scaling exponent presented in this section. 

In Fig. [HI we plot the scaling function Q + vs yi for y 2 = 8.46 at L = 180, using j e and z/ e 
to calculate Q + , and evaluating xl numerically according to Eq. ( TTTh . In addition, we plot 
in the same figure the scaling function (yi, 2/2) = X^L^ 1 ^ vs yx, for the same values of 
y 2 and L. A fit to all four Q + data points yields a scaling exponent —1.60 ± 0.03 (q = 0.02), 
while a fit to the last three Q + data points yields a scaling exponent —1.71 ±0.05 (q = 0.25). 
We will provide evidence in the next paragraph that only the last three Q + data points, 
and not the first, satisfy the asymptotic condition yi ^> y 2 . Thus, these data are consistent 
with power-law scaling of xl with exponent — 7d = — 7 C = —7/4 = —1.75. Attempts to 
fit the Qr_ data to all four and the last three data points yield nominal scaling exponents 
-1.73 ± 0.01 (q < 10~ 15 ) and -1.81 ± 0.02 (q = 2.3 x 10~ 14 ), respectively. Thus, while in 
Fig. [9] it appears that the power-law relationships with these nominal scaling exponents give 
respectable visual fits to the data, there are variations in the data which, while small, 
are larger than the statistical error bars, and which prevent a statistically significant fit. We 
will describe the causes of these variations later in this section. 

We present in Fig. [10] a plot of Q + and vs y\ at y 2 = 0, again for L = 180, for a 
larger range from y x = 30.3 to 477. With y 2 = 0, we expect that y\ = 30.3 (and indeed, 
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essentially any nonzero value of yi) should satisfy the asymptotic scaling condition y\ ^> y 2 . 
An attempt to fit to all six Q + data points in Fig. [10] gives a nominal scaling exponent 
-1.65 ± 0.03 (q = 3.8 x 10~ 7 ), while a fit to the first five points {y x = 30.3 through 280) 
gives a scaling exponent —1.74 ±0.03 (q = 0.07). Excluding the first data point at yi = 30.3 
has little effect on either fit. This supports the assumption that with y 2 = 0, the asymptotic 
scaling condition yi 3> y 2 holds for yi = 30.3, while for y 2 = 8.46, as used in Fig. [91 the 
asymptotic scaling condition does not hold for y\ = 30.3. Attempts to fit all of, and the first 
five of, the data points to power-law scaling again yield only nominal scaling exponents 
-2.01 ± 0.01 (q < 10~ 15 ) and -2.10 ± 0.01 (q < 10~ 15 ), respectively. 

We considered that the low statistical significance of the fits to the Qr: data could be 
caused by underestimation of the error bars on X®. These error bars were calculated by 
(i) finding the correlation time in the numerical data series Qi from the simulation, and 
sampling data at intervals of twice the correlation time; (ii) dividing this sampled data into 
k > 16 groups and calculating the value of X® within each group; (iii) finding the mean and 
standard error of this collection of X® values. As a check on self-consistency, we performed 
several independent calculations of X® by this method, and found that the standard error 
of these values (corrected for small sample size) was comparable to the standard error found 
within each calculation. Thus, we have strong evidence that the error bars for X® (and Q+) 
are accurate. 

These scaling results can be understood in light of the observations in Sec. IVII on the 
relationship between X® and xl- We can reasonably assume that the breakdown in scaling 
of Q + past i/i = 280 (P = 350 MCSS, 9 = 1.56) in Fig. [TOloccurs because this is the boundary 
of the critical region. The small variations of the data for P < 350 MCSS in Fig. [TPl from 
a scaling relationship with exponent — 5d = —1.75 then have three main causes. The first 
cause is the multiplication of the accurately scaling function Q + by the 9- and ^-dependent 
value T e g, according to Eq. ( fl5l) . However, such a variation would also occur in an analogous 
plot for the scaling of Xf = L 2 ((m 2 ) - (m) 2 ) vs y Xfi = eL 1 ^ = ((T - T c ) /T c ) L l / V in the 
equilibrium Ising model, since the susceptibility Xl! = d(m)/dH scales with exponent — j e , 
and Xff is related to Xl 1 by ^ ne e-dependent temperature T according to Eq. fflBT) . This 
effect is small enough to be neglected in equilibrium critical scaling, and, since the change in 
T e Q from 9 = 0.02 to 9 = 0.4 is comparable to the change in T from e = 0.02 to e = 0.4 in the 
equilibrium transition, it can also be neglected here. The second cause is the presence of the 
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nonlinear regimes in the plots of vs Xl at low xl, resulting in non-zero X^-intercepts 
in the application of Eq. ([15]) to Figs. [5] and [7J Because of this, division of the X® data 
by the appropriate T e g values (given in the caption of Fig. [5]) does not quite reproduce the 
corresponding xl data, and (even below 9 = 0.4) does not quite result in scaling consistent 
with 5^ = 1-75 with statistical signficance. The third and most signficant cause of the 
variations of the data is the 'doubly linear' behavior observed in Fig. [7J for 9 > 0.4, 
which prevents identification of a unique T e g in this range. 

The assumption that X® can be used as a proxy for xl is thus fairly well justified 
close to the critical period, where T e g varies over a limited range and the more complicated 
effects observed at 9 > 0.4 are not relevant. This is supported by Fig. [9] where the data 
points cluster closely around the line corresponding to power-law scaling with exponent 
— 1.73 ~ — 7 e . However, because of the first two causes just described, there are small 
systematic variations in the data which prevent a statistically significant fit to a pure 
power-law relationship as a function of y\. 

In order to clarify the relationship of these scaling results to those in previous work, we 

also plot in Fig.dTJdata of Q ] +\y u Jfc) = X^L~^' V = ((Q 2 ) - (\Q\) 2 ) L^'^ vs y x at y 2 = 0, 

\x\ 

again using the equilibrium values 7 e and v e to calculate y\ . This can be directly compared 



to Fig. 11(d) in Ref. 



23], in which the quantity we call X^ was called X®. Attempted fits 
to all five data points and to the last four data points of Q + in Fig. [10] produce nominal 
scaling exponents —1.60 ± 0.02 and —1.69 ± 0.02 (both with q < 10~ 15 ). The agreement in 
Fig. 11(d) of Ref. 23( of the line with slope —7/4 with the data for 9 > 9 C must therefore 
be viewed as qualitative. The method used in Ref. 2^] to numerically estimate 7d, however, 
which involves finite-size scaling at the critical period, is fully consistent with the results of 
this paper, since at each period with 9 < 0.4 we have found that the FDR in Eq. (Tl~5]) holds 
to a very good approximation. 

Finally, in Fig. [12] we plot vs y 2 at P = P c , to study its scaling in the regime y\ -C y 2 - 
As just noted, the use of X® as a proxy for xl is well justified at P = P c by our results. 
Power-law scaling is perhaps suggested in the range 20 < y 2 < 50 for L = 180, and it is 
clearly obeyed from y 2 = 8.42 to 84.2 for L = 256. The scaling exponent was determined as 
(1 — 5d)/^d = —0.914 ± 0.030, which is consistent with the corresponding equilibrium value 
(1 - 5 C )/S C = -14/15 w -0.933. 
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VIII. CONCLUSIONS AND OUTLOOK 



In this article, we have continued the computational study of the dynamic phase transition 
(DPT) in the two-dimensional kinetic Ising model exposed to a periodically oscillating field, 



which was begun in Refs 



10l . l23j . We have established two distinct but related results 
about the field conjugate to the dynamic order parameter. First, we have identified the 
period-averaged magnetic field, or 'bias field', as an important component of the full 
conjugate field. This claim is supported by numerical evidence that the dynamic order 
parameter and its susceptibility follow critical scaling with respect to H^. In particular, the 
scaling exponent 5^ of the conjugate field was determined for the first time, and found by 
finite-size scaling analysis of large-scale kinetic Monte Carlo simulations to be equal to the 
critical-isotherm exponent for the equilibrium Ising transition, 5 e = 15. Furthermore, in 
agreement with previous results {23], the dynamic scaling exponents 7d, /3d> and z/ d were also 
found to equal their equilibrium Ising counterparts, 7 e = 7/8, JL = 1/8, and — 

These results further strengthen previous numerical 0, Q, 23 1 and analytical 28, 
claims that the DPT in a periodically driven two-dimensional kinetic Ising model belongs to 
the universality class of the equilibrium two-dimensional Ising model. However, with respect 



to the direct applicability of the symmetry arguments of Refs. [32 




331 ] , we caution the reader 
that what is claimed in the present paper (as well as in Ref. 28]) is only equivalence of the 
phase transitions in the driven kinetic Ising model and the equilibrium Ising model. Outside 
the critical region, it is neither clear how closely P — P c and Hi play the roles of T — T c and 
the ordinary magnetic field, respectively, nor how closely the dynamic order parameter, Q, 
corresponds to the average equilibrium magnetization. From our discussion of the FDR in 
Sec. IVI| it appears likely that one or more of these relations break down outside the critical 
region. Much theoretical work remains to be done in this area. 

The second main result of this article is that a fluctuation-dissipation relation (FDR), 
that is, a proportionality relation between the scaled fluctuations X® = L 2 ((Q 2 ) — (Q) 2 ) 
and the susceptibility xl with a slope we have called T e ff, holds for a range of periods above 
P c and for a range of bias fields around = 0. We stress again that we have found the FDR 
of Eq. (fl5j) to hold only in the critical region in this nonequilibrium system, in contrast to the 
equilibrium FDR of Eq. ffl6l) which follows directly from the partition function, and which 
thus holds everywhere. We note that, for the parameters used in our computation at least, 
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the critical region in which the nonequilibrium FDR holds (P < 190 MCSS) is somewhat 
smaller than the critical region in which power-law scaling is obeyed (P < 350 MCSS). 
In previous work, when the conjugate field had not been identified, the scaled fluctuations 
X® were used as a proxy for the (then unknown) quantity xl- The evidence for the FDR 
presented here shows this assumption to be fully justified at the critical period (see Fig. [T21) . 
and to be a very good approximation - nearly as good as the use of the scaled fluctuations 
as a proxy for the susceptibility in the equilibrium Ising model - in the critical region where 
the FDR holds. 

There are at least three further computational projects suggested by the progress re- 
ported here. The first is to investigate whether the field functions as the conjugate field, 
with scaling exponents consistent with the equilibrium Ising transition for periods P < P c , 
below the critical period. In the equilibrium system, the study of critical scaling in nonzero 
field for T < T c is complicated by the long time correlations and strong finite-size effects 
which accompany the bimodal distributions of magnetization below T c . Similar effects would 
complicate the investigation of scaling with respect to Hi, in the DPT for P < P c , but the 
advanced techniques [4fj[ H] which make the equilibrium simulations tractable do not ex- 
tend obviously to the nonequilibrium case. The second computational project suggested is 
to determine the nature of the full conjugate field H c . The third project would be to study 
the FDR at different values of the temperature, T, and the amplitude, H , of the driving 
field. Finally, we remark that it would be very desirable to extend the current understanding 
of the theory of nonequilibrium steady states to include the conjugate field H^, the FDR 
found in the critical region, and the scaling of Hb- 
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FIG. 1: (Color online.) Log-log plot of the dynamic order parameter, (Q), vs bias field, Hb, at 
P = P c for L = 90, 128, 180, and 256. A least-squares fit to power-law scaling of the L = 256 data, 
in the range between the labels A and B above, produced a statistically significant fit with scaling 
exponent = 0.0673 ±0.0008 (corresponding to 5& = 14.85 ±0.18). The dotted line corresponds 
to the scaling exponent l/5d = 0.0673. A reference line representing scaling with the equilibrium 
Ising exponent, 5 C = 15 (l/<5 e = 0.0666), is shown as the dashed line. 
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FIG. 2: (Color online.) Log-log plots of the scaling function J~+(yi,y2) vs y\ for lattice sizes 
L = 90,128, and 180. The values of y x plotted are 4.00, 17.2, 30.3, 43.4, 69.7, 147, and (in (b)) 
280 and 477. (a) The data are shown for the values of y2 listed on the plot. The best-fit line for 
the last five points of the L = 180 data at ?/2 = 3.39, with slope —1.76 ± 0.07, is included along 
with a reference line with the slope — j e = —1.75. (b) The data for y^ = 3.39 with two additional 
yi values illustrates the boundary of the regime of power-law scaling. 
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FIG. 3: (Color online.) Log-log plot of the scaling function F+(yi,y2) y s 2/2 = J)L^I V for 
lattice sizes L = 90, 128, and 180, for the constant values of y\ labeled in the plot. The values of 
2/2 used are yi = 3.39, 8.46, 16.9, 33.9, 84.6, and 169. The dotted line represents the best fit to the 
first five points of the L = 180 data at y\ = 149, and has a slope of 1.01 ± 0.01. The dashed line 
shows the slope value of 1, expected from Eq. (Tl~3j) . 



22 




FIG. 4: (Color online.) Log-log plot of the scaling function F+(yi,y2) y s 2/2 = (Hf,/ J)L /3S ^ U for 
lattice sizes L = 128, 180, and 256, at the critical period P c , where y% = 0. In the L = 256 data, 
near the values y\ = 16, 32, 165, and 335, two closely spaced data points are actually plotted. The 
best-fit line to the L = 256 data in the range 8.46 < y% < 84.6, shown as a dotted line in the 
plot, corresponds to a scaling exponent = 0.0673 ± 0.0008. A reference line corresponding to 
scaling exponent l/5 e = 1/15 = 0.0666 is also shown. These results are in complete agreement 
with those shown in Fig. [TJ 
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FIG. 5: (Color online.) The scaled fluctuations Xj- of the dynamic order paramater plotted vs its 
susceptibility xl to the bias field Hf,, calculated at L = 180, for periods P = 140, 150, 160, 170 and 
190 MCSS. The quantity \L was calculated using the numerical derivative in Eq. (|17p . The best-fit 
lines shown, whose slopes increase monotically with the period P of the data to which they were fit, 
were calculated as (3.239 J)±l + 10.42, (3.557 J)xl +8.735, (3.980 J)±l +3.889, (4.232 J) X l + 5.868 
and (4.497J)xl + 5.371, respectively. 
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FIG. 6: (Color online.) Closeup of Fig. [5l showing the relationship of X^ and xl at low values 
of XLi which correspond to large values of the bias field H^. For P = 150, 170 and 190 MCSS, 
data have been taken (and are shown) down to very low values of xl, where the breakdown of the 
linear relationship between X® and xl can be clearly seen. The dashed lines are the same best-fit 
lines shown in Fig. [5j 
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FIG. 7: (Color online.) The scaled fluctuations of the dynamic order paramater plotted vs its 
susceptibility xl to the bias field Hb, calculated for lattice size L = 180, at periods P = 220 and 250 
MCSS. At each period, the data were fit (purely phenomenologically) to two linear relationships. 
For P = 220 MCSS, the fits were calculated as (6.265 J)xl -7.497 at low X L, and (4.161 J)xl + 19.94 
at high xl- For P = 250 MCSS, the fits were (6.485 J)xl -9.240 at low X L, and (0.2726 J)xl +67.92 
at high xl- 
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FIG. 8: (Color online.) The effective temperature T e g- , obtained as the slopes of the linear fits to 
the data in Figs. E] and plotted vs 9 = (P - P c ) /P c . For the values = 0.606 and 0.825 (P = 220 
and 250 MCSS), the slopes of both linear regimes fit in Fig. [7] are plotted as T e g values, using filled 
squares and diamonds rather than filled circles. The straight line is a weighted least-squares fit to 
the data below an 0.4 (P < 190 MCSS), and has a slope of 2. 97 J. 



27 




FIG. 9: (Color online.) Log-log plots of £/+ (2/1,2/2) and <5+ (2/1,2/2) vs 2/1, over the range 2/1 = 30.3 
through 149, for 2/2 = 8.46 at L = 180. The relatively small error bars on each data point can 
be seen inside the larger symbols. The solid and dash-dash-dotted lines are fits to all four and 
the last three Q+ data points, respectively, and correspond to scaling exponents —1.60 ± 0.03 and 
— 1.71±0.05. The dotted and dash-dotted lines are the result of attempts to fit all four and the last 
three Q¥ data points, respectively. They correspond to nominal scaling exponents —1.73 ± 0.01 
and —1.81 ± 0.02. The dashed line is a reference line corresponding to scaling exponent —1.75. 
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FIG. 10: (Color online.) Log-log plots of G+(yi,y2) and G+(yi,U2) vs yx, over the range y\ = 30.3 
through 477, for y% = at lattice size L = 180. The relatively small error bars on each data 
point can be seen inside the larger symbols. The solid and dash-dash-dotted lines show the result 
of attempts to fit all six and the first five G+ data points, and correspond to a nominal scaling 
exponent —1.65 ± 0.03, and a statistically significant scaling exponent —1.74 ± 0.03, respectively. 
The dotted and dash-dotted lines are the result of attempts to fit all six and the first five data 
points, and correspond to nominal scaling exponents —2.01 ± 0.01 and —2.10 ± 0.01, respectively. 
The dashed line is a reference line corresponding to scaling exponent —1.75. 
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FIG. 11: (Color online.) Log-log plot of Q + (2/1 5 2/2 ) v s Ui, over the range y\ = 30.3 through 
149, for y2 = at lattice size L = 180. The relatively small error bars on each data point can 
be seen inside the larger symbols. The solid and dotted lines are the results of attempts to fit 
all five and the last four data points with a power-law relationship, and correspond to nominal 
scaling exponents of —1.59 ± 0.02 and —1.70 ± 0.03, respectively. The dashed line is a reference 
line corresponding to a scaling exponent of —1.75. 
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FIG. 12: (Color online.) Log-log plot of g£(|ft, jfe) vs y 2 = {H b /J)L^ 5 / U for lattice sizes L = 128, 
180, and 256, at the critical period P c , where y\ = 0. The best-fit line to the L = 256 data in 
the range 8.46 < yi < 84.6, shown as a dotted line in the plot, corresponds to a scaling exponent 
(1 — 5d)/5d = —0.914 ± 0.029. A reference line (dashed), corresponding to a scaling exponent 
(1 - 5 e )/5 e = -14/15 w -0.933, is also shown. 
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